//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : IncomingWaves.cpp
//----------------------------------------------------------------------------

#include <math.h>
#ifndef _WINDOWS
#include "mermaid.h"
#endif
#include "FloatUtils.h"
#include "Array1D.h"
#include "Array2D.h"
#include "IncomingWaves.h"
#include "assertplus.h"

#ifdef _WINDOWS
#pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data
#pragma warning(disable: 4305) // truncation from 'double' to 'float'
#endif

#define NUMBER   wvin_spec_count + 1
#define COLUMNS  (unsigned)(Nx1 + 3)
#define ROWS     (unsigned)(Nx2 + 3)

/// constructor
IncomingWaves::IncomingWaves(int Nx1,
                             int Nx2,
                             REAL dx1,
                             REAL dx2,
                             const Array2D<REAL> &h0,
                             const Array2D<REAL> &U1,
                             const Array2D<REAL> &U2,
                             const Array2D<SHORT> &wet,
                             const Array2D<REAL> &kappa,
                             const Array2D<REAL> &vD0,
                             const Array2D<REAL> &vM0,
                             const Array2D<REAL> &vN0,
                             int wvin_spec_count,
                             const Array1D<REAL> &wvin_spec_ampl,
                             const Array1D<REAL> &wvin_spec_phase,
                             const Array1D<REAL> &wvin_spec_freq,
                             const Array1D<REAL> &wvin_spec_dir,
                             REAL g,
                             VariationalModel varmod) 
: m_Nx1(Nx1)
, m_Nx2(Nx2)
, m_dx1(dx1)
, m_dx2(dx2)
, m_g(g)
, m_varmod(varmod)
, m_h0(h0)
, m_U1(U1)
, m_U2(U2)
, m_wet(wet)
, m_kappa(kappa)
, m_vD0(vD0)
, m_vN0(vN0)
, m_wvin_spec_count(wvin_spec_count)
, m_wvin_spec_ampl(wvin_spec_ampl)
, m_wvin_spec_phase(wvin_spec_phase)
, m_wvin_spec_freq(wvin_spec_freq)
, m_wvin_spec_dir(wvin_spec_dir)
, m_wvin_N_count(COLUMNS)
, m_wvin_S_count(COLUMNS)
, m_wvin_E_count(ROWS)
, m_wvin_W_count(ROWS)
, m_wvin_N_cen_phase(NUMBER, COLUMNS)
, m_wvin_S_cen_phase(NUMBER, COLUMNS)
, m_wvin_N_int_phase(NUMBER, COLUMNS)
, m_wvin_S_int_phase(NUMBER, COLUMNS)
, m_wvin_N_ext_phase(NUMBER, COLUMNS)
, m_wvin_S_ext_phase(NUMBER, COLUMNS)
, m_wvin_N_cen_facphi(NUMBER, COLUMNS)
, m_wvin_S_cen_facphi(NUMBER, COLUMNS)
, m_wvin_N_cen_facpsi(NUMBER, COLUMNS)
, m_wvin_S_cen_facpsi(NUMBER, COLUMNS)
, m_wvin_N_int_facphi(NUMBER, COLUMNS)
, m_wvin_S_int_facphi(NUMBER, COLUMNS)
, m_wvin_N_ext_facphi(NUMBER, COLUMNS)
, m_wvin_S_ext_facphi(NUMBER, COLUMNS)
, m_wvin_N_int_facpsi(NUMBER, COLUMNS)
, m_wvin_S_int_facpsi(NUMBER, COLUMNS)
, m_wvin_N_ext_facpsi(NUMBER, COLUMNS)
, m_wvin_S_ext_facpsi(NUMBER, COLUMNS)
, m_wvin_N_cen_ampl(NUMBER, COLUMNS)
, m_wvin_N_cen_freq(NUMBER, COLUMNS)
, m_wvin_S_cen_ampl(NUMBER, COLUMNS)
, m_wvin_S_cen_freq(NUMBER, COLUMNS)
, m_wvin_N_int_ampl(NUMBER, COLUMNS)
, m_wvin_N_int_freq(NUMBER, COLUMNS)
, m_wvin_S_int_ampl(NUMBER, COLUMNS)
, m_wvin_N_ext_ampl(NUMBER, COLUMNS)
, m_wvin_N_ext_freq(NUMBER, COLUMNS)
, m_wvin_S_ext_ampl(NUMBER, COLUMNS)
, m_wvin_S_int_freq(NUMBER, COLUMNS)
, m_wvin_S_ext_freq(NUMBER, COLUMNS)
, m_wvin_E_cen_phase(NUMBER, ROWS)
, m_wvin_W_cen_phase(NUMBER, ROWS)
, m_wvin_E_int_phase(NUMBER, ROWS)
, m_wvin_E_ext_phase(NUMBER, ROWS)
, m_wvin_W_int_phase(NUMBER, ROWS)
, m_wvin_W_ext_phase(NUMBER, ROWS)
, m_wvin_E_cen_facphi(NUMBER, ROWS)
, m_wvin_E_cen_facpsi(NUMBER, ROWS)
, m_wvin_W_cen_facphi(NUMBER, ROWS)
, m_wvin_E_int_facphi(NUMBER, ROWS)
, m_wvin_E_ext_facphi(NUMBER, ROWS)
, m_wvin_W_cen_facpsi(NUMBER, ROWS)
, m_wvin_E_int_facpsi(NUMBER, ROWS)
, m_wvin_W_int_facphi(NUMBER, ROWS)
, m_wvin_E_ext_facpsi(NUMBER, ROWS)
, m_wvin_W_ext_facphi(NUMBER, ROWS)
, m_wvin_W_int_facpsi(NUMBER, ROWS)
, m_wvin_W_ext_facpsi(NUMBER, ROWS)
, m_wvin_E_cen_ampl(NUMBER, ROWS)
, m_wvin_E_cen_freq(NUMBER, ROWS)
, m_wvin_W_cen_ampl(NUMBER, ROWS)
, m_wvin_E_int_ampl(NUMBER, ROWS)
, m_wvin_W_cen_freq(NUMBER, ROWS)
, m_wvin_E_int_freq(NUMBER, ROWS)
, m_wvin_E_ext_ampl(NUMBER, ROWS)
, m_wvin_E_ext_freq(NUMBER, ROWS)
, m_wvin_W_int_ampl(NUMBER, ROWS)
, m_wvin_W_int_freq(NUMBER, ROWS)
, m_wvin_W_ext_ampl(NUMBER, ROWS)
, m_wvin_W_ext_freq(NUMBER, ROWS)
#ifdef USE_CUDAINCOMINGWAVES
, m_pCudaIncomingWaves(0) // Class CudaIncomingWaves
#endif
{
#ifdef USE_CUDAINCOMINGWAVES
    m_pCudaIncomingWaves = new CudaIncomingWaves(m_Nx1, m_Nx2, m_wvin_spec_count);
#endif

    initialize(vM0);

#ifdef USE_CUDAINCOMINGWAVES
    m_pCudaIncomingWaves->initialize(&m_wvin_N_count,
                                     &m_wvin_N_ext_ampl,
                                     &m_wvin_N_ext_phase,
                                     &m_wvin_N_ext_freq,
                                     &m_wvin_N_ext_facphi,
                                     &m_wvin_N_ext_facpsi,
                                     &m_wvin_N_cen_ampl,
                                     &m_wvin_N_cen_phase,
                                     &m_wvin_N_cen_freq,
                                     &m_wvin_N_cen_facphi,
                                     &m_wvin_N_cen_facpsi,
                                     &m_wvin_N_int_ampl,
                                     &m_wvin_N_int_phase,
                                     &m_wvin_N_int_freq,
                                     &m_wvin_N_int_facphi,
                                     &m_wvin_N_int_facpsi,
                                     &m_wvin_E_count,
                                     &m_wvin_E_ext_ampl,
                                     &m_wvin_E_ext_phase,
                                     &m_wvin_E_ext_freq,
                                     &m_wvin_E_ext_facphi,
                                     &m_wvin_E_ext_facpsi,
                                     &m_wvin_E_cen_ampl,
                                     &m_wvin_E_cen_phase,
                                     &m_wvin_E_cen_freq,
                                     &m_wvin_E_cen_facphi,
                                     &m_wvin_E_cen_facpsi,
                                     &m_wvin_E_int_ampl,
                                     &m_wvin_E_int_phase,
                                     &m_wvin_E_int_freq,
                                     &m_wvin_E_int_facphi,
                                     &m_wvin_E_int_facpsi,
                                     &m_wvin_S_count,
                                     &m_wvin_S_ext_ampl,
                                     &m_wvin_S_ext_phase,
                                     &m_wvin_S_ext_freq,
                                     &m_wvin_S_ext_facphi,
                                     &m_wvin_S_ext_facpsi,
                                     &m_wvin_S_cen_ampl,
                                     &m_wvin_S_cen_phase,
                                     &m_wvin_S_cen_freq,
                                     &m_wvin_S_cen_facphi,
                                     &m_wvin_S_cen_facpsi,
                                     &m_wvin_S_int_ampl,
                                     &m_wvin_S_int_phase,
                                     &m_wvin_S_int_freq,
                                     &m_wvin_S_int_facphi,
                                     &m_wvin_S_int_facpsi,
                                     &m_wvin_W_count,
                                     &m_wvin_W_ext_ampl,
                                     &m_wvin_W_ext_phase,
                                     &m_wvin_W_ext_freq,
                                     &m_wvin_W_ext_facphi,
                                     &m_wvin_W_ext_facpsi,
                                     &m_wvin_W_cen_ampl,
                                     &m_wvin_W_cen_phase,
                                     &m_wvin_W_cen_freq,
                                     &m_wvin_W_cen_facphi,
                                     &m_wvin_W_cen_facpsi,
                                     &m_wvin_W_int_ampl,
                                     &m_wvin_W_int_phase,
                                     &m_wvin_W_int_freq,
                                     &m_wvin_W_int_facphi,
                                     &m_wvin_W_int_facpsi);
#endif
}

/// destructor
IncomingWaves::~IncomingWaves()
{
}

/// initializes boundary conditions
void IncomingWaves::initialize(const Array2D<REAL> &vM0)
{
    REAL pi = 4.0 * atan(1.0);

    // North boundary 
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        m_wvin_N_count[i1] = 0;
        for (int wv = 1; wv <= m_wvin_spec_count; ++wv)
        {
            if ((m_wvin_spec_dir[wv] >= 1.0 * pi && m_wvin_spec_dir[wv] <= 2.0 * pi) ||
                (m_wvin_spec_dir[wv] == 0.0 * pi))
            {
                ++m_wvin_N_count[i1];
                int j_now = m_wvin_N_count[i1];
                // the exterior points
                int i2 = m_Nx2 + 1;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_N_ext_ampl[j_now][i1]), &(m_wvin_N_ext_phase[j_now][i1]),
                          &(m_wvin_N_ext_freq[j_now][i1]), &(m_wvin_N_ext_facphi[j_now][i1]),
                          &(m_wvin_N_ext_facpsi[j_now][i1]));
                // the central points
                i2 = m_Nx2;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_N_cen_ampl[j_now][i1]), &(m_wvin_N_cen_phase[j_now][i1]),
                          &(m_wvin_N_cen_freq[j_now][i1]), &(m_wvin_N_cen_facphi[j_now][i1]),
                          &(m_wvin_N_cen_facpsi[j_now][i1]));
                // the interior points
                i2 = m_Nx2 - 1;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_N_int_ampl[j_now][i1]), &(m_wvin_N_int_phase[j_now][i1]),
                          &(m_wvin_N_int_freq[j_now][i1]), &(m_wvin_N_int_facphi[j_now][i1]),
                          &(m_wvin_N_int_facpsi[j_now][i1]));
            }
        }
    }

    // East boundary 
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        m_wvin_E_count[i2] = 0;
        for (int wv = 1; wv <= m_wvin_spec_count; ++wv)
        {
            if (m_wvin_spec_dir[wv] >= 0.5 * pi && m_wvin_spec_dir[wv] <= 1.5 * pi)
            {
                ++m_wvin_E_count[i2];
                int j_now = m_wvin_E_count[i2];
                // the exterior points
                int i1 = m_Nx1 + 1;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_E_ext_ampl[j_now][i2]), &(m_wvin_E_ext_phase[j_now][i2]),
                          &(m_wvin_E_ext_freq[j_now][i2]), &(m_wvin_E_ext_facphi[j_now][i2]),
                          &(m_wvin_E_ext_facpsi[j_now][i2]));
                // the central points
                i1 = m_Nx1;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_E_cen_ampl[j_now][i2]), &(m_wvin_E_cen_phase[j_now][i2]),
                          &(m_wvin_E_cen_freq[j_now][i2]), &(m_wvin_E_cen_facphi[j_now][i2]),
                          &(m_wvin_E_cen_facpsi[j_now][i2]));
                // the interior points
                i1 = m_Nx1 - 1;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_E_int_ampl[j_now][i2]), &(m_wvin_E_int_phase[j_now][i2]),
                          &(m_wvin_E_int_freq[j_now][i2]), &(m_wvin_E_int_facphi[j_now][i2]),
                          &(m_wvin_E_int_facpsi[j_now][i2]));
            }
        }
    }

    // South boundary 
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        m_wvin_S_count[i1] = 0;
        for (int wv = 1; wv <= m_wvin_spec_count; ++wv)
        {
            if ((m_wvin_spec_dir[wv] >= 0.0 * pi && m_wvin_spec_dir[wv] <= pi) ||
                m_wvin_spec_dir[wv] == 2.0 * pi)
            {
                ++m_wvin_S_count[i1];
                int j_now = m_wvin_S_count[i1];
                // the exterior points
                int i2 = 0;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_S_ext_ampl[j_now][i1]), &(m_wvin_S_ext_phase[j_now][i1]),
                          &(m_wvin_S_ext_freq[j_now][i1]), &(m_wvin_S_ext_facphi[j_now][i1]),
                          &(m_wvin_S_ext_facpsi[j_now][i1]));
                // the central points
                i2 = 1;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_S_cen_ampl[j_now][i1]), &(m_wvin_S_cen_phase[j_now][i1]),
                          &(m_wvin_S_cen_freq[j_now][i1]), &(m_wvin_S_cen_facphi[j_now][i1]),
                          &(m_wvin_S_cen_facpsi[j_now][i1]));
                // the interior points
                i2 = 2;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_S_int_ampl[j_now][i1]), &(m_wvin_S_int_phase[j_now][i1]),
                          &(m_wvin_S_int_freq[j_now][i1]), &(m_wvin_S_int_facphi[j_now][i1]),
                          &(m_wvin_S_int_facpsi[j_now][i1]));
            }
        }
    }

    // West boundary 
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        m_wvin_W_count[i2] = 0;
        for (int wv = 1; wv <= m_wvin_spec_count; ++wv)
        {
            if ((m_wvin_spec_dir[wv] >= 0.0 * pi && m_wvin_spec_dir[wv] <= 0.5 * pi) ||
                (m_wvin_spec_dir[wv] >= 1.5 * pi && m_wvin_spec_dir[wv] <= 2.0 * pi))
            {
                ++m_wvin_W_count[i2];
                int j_now = m_wvin_W_count[i2];
                // the exterior points 
                int i1 = 0;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_W_ext_ampl[j_now][i2]), &(m_wvin_W_ext_phase[j_now][i2]),
                          &(m_wvin_W_ext_freq[j_now][i2]), &(m_wvin_W_ext_facphi[j_now][i2]),
                          &(m_wvin_W_ext_facpsi[j_now][i2]));
                // the central points 
                i1 = 1;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_W_cen_ampl[j_now][i2]), &(m_wvin_W_cen_phase[j_now][i2]),
                          &(m_wvin_W_cen_freq[j_now][i2]), &(m_wvin_W_cen_facphi[j_now][i2]),
                          &(m_wvin_W_cen_facpsi[j_now][i2]));
                // the interior points 
                i1 = 2;
                wvin_init(i1, i2, wv, vM0, &(m_wvin_W_int_ampl[j_now][i2]), &(m_wvin_W_int_phase[j_now][i2]),
                          &(m_wvin_W_int_freq[j_now][i2]), &(m_wvin_W_int_facphi[j_now][i2]),
                          &(m_wvin_W_int_facpsi[j_now][i2]));
            }
        }
    }
}

void IncomingWaves::startup(REAL t, bool maxDepth, Array2D<REAL>* a_pZeta, Array2D<REAL>* a_pPhi)
{
    *a_pZeta = 0.0;
    *a_pPhi = 0.0;
    REAL hmax = 0.0;
    REAL kappa = 0.0;
    if (maxDepth)
    {
        // find maximum depth with corresponding kappa
        for (int i1 = 0; i1 < (int)m_h0.dim1(); ++i1)
        {
            for (int i2 = 0; i2 < (int)m_h0.dim2(); ++i2)
            {
                if (m_h0[i1][i2] > hmax)
                {
                    hmax = m_h0[i1][i2];
                    kappa = m_kappa[i1][i2];
                }
            }
        }
    }
    for (int wv = 1; wv <= m_wvin_spec_count; wv++)
    {
        REAL cosine = cos(m_wvin_spec_dir[wv]);
        REAL sine = sin(m_wvin_spec_dir[wv]);
        for (int i1 = 0; i1 < (int)a_pZeta->dim1(); ++i1)
        {
            for (int i2 = 0; i2 < (int)a_pZeta->dim2(); ++i2)
            {
                REAL U = m_U1[i1][i2] * cosine + m_U2[i1][i2] * sine;
                REAL k;
                if (maxDepth) // maximum depth and constant kappa
                {
                    vbm_disper(&k, m_wvin_spec_freq[wv], hmax, U, kappa);
                }
                else // depth and kappa per grid point
                {
                    vbm_disper(&k, m_wvin_spec_freq[wv], m_h0[i1][i2], U, m_kappa[i1][i2]);
                }
                REAL k1 = k * cosine;
                REAL k2 = k * sine;
                //?                REAL sigma = m_wvin_spec_freq[wv] - k * U;
                if (m_h0[i1][i2] > 0.0 && m_wet[i1][i2])
                {
                    REAL x1 = (i1 - 1) * m_dx1;
                    REAL x2 = (i2 - 1) * m_dx2;
                    REAL phase = k1 * x1 + k2 * x2 - m_wvin_spec_phase[wv];
                    REAL omega = m_wvin_spec_freq[wv];
                    REAL ampl = m_wvin_spec_ampl[wv];
                    REAL factor = omega * t - phase;
                    (*a_pZeta)[i1][i2] += ampl * sin(factor);
                    int fudge = 1;
                    if (i1 == 0 || i1 == (int)a_pZeta->dim1() - 1 ||
                        i2 == 0 || i2 == (int)a_pZeta->dim2() - 1)
                    {
                        fudge = 9;
                    }
                    else
                    {
                        for (int j1 = i1 - 1; j1 <= i1 + 1; ++j1)
                        {
                            for (int j2 = i2 - 1; j2 <= i2 + 1; ++j2)
                            {
                                if (!m_wet[j1][j2])
                                {
                                    ++fudge;
                                }
                            }
                        }
                    }
                    (*a_pPhi)[i1][i2] += ampl * omega / k * cos(factor) / fudge;
                }
            }
        }
    }
    //    PrintArray2Dcontents("",(*a_pZeta),0,a_pZeta->dim1()-1,0,a_pZeta->dim2()-1);
    //    PrintArray2Dcontents("",(*a_pPhi),0,a_pZeta->dim1()-1,0,a_pZeta->dim2()-1);
}



/// computes boundary parameters related to the incoming waves
void IncomingWaves::wvin_init(int i1,
                              int i2,
                              int wv,
                              const Array2D<REAL> &vM0,
                              REAL *pwvin_W_ext_ampl,
                              REAL *pwvin_W_ext_phase,
                              REAL *pwvin_W_ext_freq,
                              REAL *pwvin_W_ext_facphi,
                              REAL *pwvin_W_ext_facpsi)
{
    REAL cosine = cos(m_wvin_spec_dir[wv]);
    REAL sine = sin(m_wvin_spec_dir[wv]);

    // compute the wave number 
    REAL U = m_U1[i1][i2] * cosine + m_U2[i1][i2] * sine;
    REAL k;
    vbm_disper(&k, m_wvin_spec_freq[wv], m_h0[i1][i2], U, m_kappa[i1][i2]);
    REAL k1 = k * cosine;
    REAL k2 = k * sine;
    REAL sigma = m_wvin_spec_freq[wv] - k * U;
    if (m_h0[i1][i2] > 0.0 && m_wet[i1][i2])
    {
        REAL x1 = (i1 - 1) * m_dx1;
        REAL x2 = (i2 - 1) * m_dx2;
        *pwvin_W_ext_phase = k1 * x1 + k2 * x2 - m_wvin_spec_phase[wv];
        *pwvin_W_ext_ampl = m_wvin_spec_ampl[wv];
        *pwvin_W_ext_freq = m_wvin_spec_freq[wv];
        *pwvin_W_ext_facphi = m_g / sigma;
        *pwvin_W_ext_facpsi = m_g / sigma * FloatUtils::square(k) *
            m_h0[i1][i2] * m_vD0[i1][i2] /
            (vM0[i1][i2] + FloatUtils::square(k) * m_vN0[i1][i2]);
    }
    else
    {
        *pwvin_W_ext_phase = 0.0;
        *pwvin_W_ext_ampl = 0.0;
        *pwvin_W_ext_freq = 1.0;
        *pwvin_W_ext_facphi = 0.0;
        *pwvin_W_ext_facpsi = 0.0;
    }
}

/// computes boundary conditions related to the incoming waves
void IncomingWaves::compute(REAL t, REAL pslow,
                            Array1D<REAL> &zeta_in_N_ext, Array1D<REAL> &phi_in_N_ext, Array1D<REAL> &psi_in_N_ext,
                            Array1D<REAL> &zeta_in_N_cen, Array1D<REAL> &phi_in_N_cen, Array1D<REAL> &psi_in_N_cen,
                            Array1D<REAL> &zeta_in_N_int, Array1D<REAL> &phi_in_N_int, Array1D<REAL> &psi_in_N_int,
                            Array1D<REAL> &zeta_in_E_ext, Array1D<REAL> &phi_in_E_ext, Array1D<REAL> &psi_in_E_ext,
                            Array1D<REAL> &zeta_in_E_cen, Array1D<REAL> &phi_in_E_cen, Array1D<REAL> &psi_in_E_cen,
                            Array1D<REAL> &zeta_in_E_int, Array1D<REAL> &phi_in_E_int, Array1D<REAL> &psi_in_E_int,
                            Array1D<REAL> &zeta_in_S_ext, Array1D<REAL> &phi_in_S_ext, Array1D<REAL> &psi_in_S_ext,
                            Array1D<REAL> &zeta_in_S_cen, Array1D<REAL> &phi_in_S_cen, Array1D<REAL> &psi_in_S_cen,
                            Array1D<REAL> &zeta_in_S_int, Array1D<REAL> &phi_in_S_int, Array1D<REAL> &psi_in_S_int,
                            Array1D<REAL> &zeta_in_W_ext, Array1D<REAL> &phi_in_W_ext, Array1D<REAL> &psi_in_W_ext,
                            Array1D<REAL> &zeta_in_W_cen, Array1D<REAL> &phi_in_W_cen, Array1D<REAL> &psi_in_W_cen,
                            Array1D<REAL> &zeta_in_W_int, Array1D<REAL> &phi_in_W_int, Array1D<REAL> &psi_in_W_int)
{
#ifdef USE_CUDAINCOMINGWAVES
    m_pCudaIncomingWaves->wvin_comp(t, pslow, 
                                    &zeta_in_N_ext, &phi_in_N_ext, &psi_in_N_ext,
                                    &zeta_in_N_cen, &phi_in_N_cen, &psi_in_N_cen,
                                    &zeta_in_N_int, &phi_in_N_int, &psi_in_N_int,
                                    &zeta_in_E_ext, &phi_in_E_ext, &psi_in_E_ext,
                                    &zeta_in_E_cen, &phi_in_E_cen, &psi_in_E_cen,
                                    &zeta_in_E_int, &phi_in_E_int, &psi_in_E_int,
                                    &zeta_in_S_ext, &phi_in_S_ext, &psi_in_S_ext,
                                    &zeta_in_S_cen, &phi_in_S_cen, &psi_in_S_cen,
                                    &zeta_in_S_int, &phi_in_S_int, &psi_in_S_int,
                                    &zeta_in_W_ext, &phi_in_W_ext, &psi_in_W_ext,
                                    &zeta_in_W_cen, &phi_in_W_cen, &psi_in_W_cen,
                                    &zeta_in_W_int, &phi_in_W_int, &psi_in_W_int);
#else
    // this is done inside the WavesComputer::compute() parallelization:
    // compute the incoming waves (4 parallel sections)

    // three function calls per section
    { //====================================================================
        { //----------------------------------------------------------------
            // north: exterior, central and interior points near boundary 
            wvin_comp(m_Nx1, m_wvin_N_count, t, pslow,
                      zeta_in_N_ext, phi_in_N_ext, psi_in_N_ext,
                      m_wvin_N_ext_ampl, m_wvin_N_ext_phase, m_wvin_N_ext_freq,
                      m_wvin_N_ext_facphi, m_wvin_N_ext_facpsi);
            wvin_comp(m_Nx1, m_wvin_N_count, t, pslow,
                      zeta_in_N_cen, phi_in_N_cen, psi_in_N_cen,
                      m_wvin_N_cen_ampl, m_wvin_N_cen_phase, m_wvin_N_cen_freq,
                      m_wvin_N_cen_facphi, m_wvin_N_cen_facpsi);
            wvin_comp(m_Nx1, m_wvin_N_count, t, pslow,
                      zeta_in_N_int, phi_in_N_int, psi_in_N_int,
                      m_wvin_N_int_ampl, m_wvin_N_int_phase, m_wvin_N_int_freq,
                      m_wvin_N_int_facphi, m_wvin_N_int_facpsi);
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            // east: exterior, central and interior points near boundary 
            wvin_comp(m_Nx2, m_wvin_E_count, t, pslow,
                      zeta_in_E_ext, phi_in_E_ext, psi_in_E_ext,
                      m_wvin_E_ext_ampl, m_wvin_E_ext_phase, m_wvin_E_ext_freq,
                      m_wvin_E_ext_facphi, m_wvin_E_ext_facpsi);
            wvin_comp(m_Nx2, m_wvin_E_count, t, pslow,
                      zeta_in_E_cen, phi_in_E_cen, psi_in_E_cen,
                      m_wvin_E_cen_ampl, m_wvin_E_cen_phase, m_wvin_E_cen_freq,
                      m_wvin_E_cen_facphi, m_wvin_E_cen_facpsi);
            wvin_comp(m_Nx2, m_wvin_E_count, t, pslow,
                      zeta_in_E_int, phi_in_E_int, psi_in_E_int,
                      m_wvin_E_int_ampl, m_wvin_E_int_phase, m_wvin_E_int_freq,
                      m_wvin_E_int_facphi, m_wvin_E_int_facpsi);
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            // south: exterior, central and interior points near boundary 
            wvin_comp(m_Nx1, m_wvin_S_count, t, pslow,
                      zeta_in_S_ext, phi_in_S_ext, psi_in_S_ext,
                      m_wvin_S_ext_ampl, m_wvin_S_ext_phase, m_wvin_S_ext_freq,
                      m_wvin_S_ext_facphi, m_wvin_S_ext_facpsi);
            wvin_comp(m_Nx1, m_wvin_S_count, t, pslow,
                      zeta_in_S_cen, phi_in_S_cen, psi_in_S_cen,
                      m_wvin_S_cen_ampl, m_wvin_S_cen_phase, m_wvin_S_cen_freq,
                      m_wvin_S_cen_facphi, m_wvin_S_cen_facpsi);
            wvin_comp(m_Nx1, m_wvin_S_count, t, pslow,
                      zeta_in_S_int, phi_in_S_int, psi_in_S_int,
                      m_wvin_S_int_ampl, m_wvin_S_int_phase, m_wvin_S_int_freq,
                      m_wvin_S_int_facphi, m_wvin_S_int_facpsi);
        } //section---------------------------------------------------------
        { //----------------------------------------------------------------
            // west: exterior, central and interior points near boundary 
            wvin_comp(m_Nx2, m_wvin_W_count, t, pslow,
                      zeta_in_W_ext, phi_in_W_ext, psi_in_W_ext,
                      m_wvin_W_ext_ampl, m_wvin_W_ext_phase, m_wvin_W_ext_freq,
                      m_wvin_W_ext_facphi, m_wvin_W_ext_facpsi);
            wvin_comp(m_Nx2, m_wvin_W_count, t, pslow,
                      zeta_in_W_cen, phi_in_W_cen, psi_in_W_cen,
                      m_wvin_W_cen_ampl, m_wvin_W_cen_phase, m_wvin_W_cen_freq,
                      m_wvin_W_cen_facphi, m_wvin_W_cen_facpsi);
            wvin_comp(m_Nx2, m_wvin_W_count, t, pslow,
                      zeta_in_W_int, phi_in_W_int, psi_in_W_int,
                      m_wvin_W_int_ampl, m_wvin_W_int_phase, m_wvin_W_int_freq,
                      m_wvin_W_int_facphi, m_wvin_W_int_facpsi);
        } //section---------------------------------------------------------

    } //====================================================================

#endif // USE_CUDAINCOMINGWAVES
}

/// computes boundary conditions for one direction and zone (ext,cen,int) of incoming waves
void IncomingWaves::wvin_comp(int Nx2,
                              const Array1D<int> &wvin_W_count,
                              REAL t,
                              REAL pslow,
                              Array1D<REAL> &zeta_in_W_ext,
                              Array1D<REAL> &phi_in_W_ext,
                              Array1D<REAL> &psi_in_W_ext,
                              const Array2D<REAL> &wvin_W_ext_ampl,
                              const Array2D<REAL> &wvin_W_ext_phase,
                              const Array2D<REAL> &wvin_W_ext_freq,
                              const Array2D<REAL> &wvin_W_ext_facphi,
                              const Array2D<REAL> &wvin_W_ext_facpsi)
{
    // this is done inside the WavesComputer::compute() parallelization:
    // called twelve times from IncomingWaves::compute()
    // compute zeta, phi and psi per wind direction (parallel for)

    // (Nx2 + 2) * wvin_W_count[i2] calculations
    for (int i2 = 0; i2 <= Nx2 + 1; ++i2)
    {
        // save some cycles by remembering pointers
        REAL* pzeta = &(zeta_in_W_ext[i2]);
        REAL* pphi  = &(phi_in_W_ext[i2]);
        REAL* ppsi  = &(psi_in_W_ext[i2]);
        *pzeta = 0.0;
        *pphi  = 0.0;
        *ppsi  = 0.0;
        int count = wvin_W_count[i2];
        for (int wv = 1; wv <= count; ++wv)
        {
            REAL phase  = wvin_W_ext_phase[wv][i2] - wvin_W_ext_freq[wv][i2] * t;
            REAL ampl   = wvin_W_ext_ampl[wv][i2];
//#define SINCOS
#ifdef SINCOS
            // proved to be NO faster!
            double cosphase, sinphase;
            sincos(phase, &sinphase, &cosphase);
            REAL cosamp = ampl * cosphase;
            REAL sinamp = ampl * sinphase;
#else            
            REAL cosamp = ampl * cos(phase);
            REAL sinamp = ampl * sin(phase);
#endif            
            *pzeta += cosamp;
            *pphi  += wvin_W_ext_facphi[wv][i2] * sinamp;
            *ppsi  += wvin_W_ext_facpsi[wv][i2] * sinamp;
        }
        *pzeta *= pslow;
        *pphi  *= pslow;
        *ppsi  *= pslow;
    }
    // for--------------------------------------------------------------
}

/// computes wave number for the linearised variational Boussinesq model
void IncomingWaves::vbm_disper(REAL *pk,
                               REAL omega,
                               REAL h,
                               REAL U,
                               REAL kappa)
{
    const REAL tol = 1.0e-5;
    const int maxIter = 20; // maximum number of iterations in Newton process
#ifdef nodig
    const int solMeth = 3; // zero-finding equation to be used in Newton process
#endif

    // compute non-dimensional Ome, Fr and kappa_h 
    REAL Ome = omega * sqrt(h / m_g);
    REAL Fr = U / sqrt(m_g * h);
    REAL kappa_h = kappa * h;

    // initialise Newton process parameters 

    // compute the parameters from the VBM 
    REAL S = 0.0;
    REAL C = 0.0;
    REAL D = 0.0;
    REAL M = 0.0;
    REAL N = 0.0;
    switch (m_varmod)
    {
        case coshypVM: // hyperbolic-cosine structure model
            S = sinh(kappa_h);
            C = cosh(kappa_h);
            D = C - S / kappa_h;
            M = 0.5 * kappa_h * (S * C - kappa_h);
            N = 0.5 / kappa_h * (2 * kappa_h * FloatUtils::square(C) + kappa_h - 3 * S * C);
            break;
        case parabolicVM: // parabolic structure model
            // NOTE ASR: this leads to very strange results, in the original Fortran code too!
            S = 0.0;
            C = 0.0;
            D = 1.0 / 3.0;
            M = 1.0 / 3.0;
            N = 2.0 / 15.0;
            break;
        default:
            ASSERTFALSE(); // Unknown vertical structure model
    }

    // initial guess for the relative wave number k*h 
    REAL kh = Ome;
    REAL Dsquare = FloatUtils::square(D);
    if (fabs(Ome) < 1.0)
    {
        if (Fr != -1.0)
        {
            kh = Ome / (Fr + 1.0);
        }
    }
    else
    {
        if (Fr != -sqrt(1.0 - Dsquare / N))
        {
            kh = Ome / (Fr + sqrt(1.0 - Dsquare / N));
        }
    }

    // Newton iteration process 
    REAL err = fabs(Ome);
    int iter = 0;
    while (fabs(err) > tol * fabs(Ome) && iter < maxIter)
    {
        ++iter;
        // f is the function to find the root from, and df is its derivative 
        REAL khsquare = FloatUtils::square(kh);
        REAL f, df;
#ifdef nodig
        if (solMeth == 1)
        {
            f = FloatUtils::square(Ome - kh * Fr) * (M + khsquare * N) -
                khsquare * (M + khsquare * (N - Dsquare));
            df = 2 * kh * (FloatUtils::square(Fr) - 1) * (M + 2 * khsquare * N)
                + 4 * FloatUtils::third(kh) * Dsquare
                - 2 * Fr * (M + 3 * khsquare * N) * Ome + 2 * kh * N * FloatUtils::square(Ome);
        }
        else if (solMeth == 2)
        {
            f = FloatUtils::square(Ome - kh * Fr) - khsquare * (M + khsquare *
                                                                (N - Dsquare)) / (M + khsquare * N);
            df = -2 * Ome * Fr + 2 * kh * (FloatUtils::square(Fr) - 1) +
                2 * FloatUtils::third(kh) * Dsquare *
                (2 * M + khsquare * N) / FloatUtils::square(M + khsquare * N);
        }
        else if (solMeth == 3)
        {
#endif // nodig
            REAL tmp1 = sqrt((M + khsquare * (N - Dsquare)) /
                             (M + khsquare * N));
            f = Ome - kh * Fr - kh * tmp1;
            df = -Fr - tmp1 + khsquare * M *
                Dsquare / tmp1 / FloatUtils::square(M + khsquare * N);
#ifdef nodig
        }
        else
        {
            ASSERT(false); // Wrong solution number
        }
#endif // nodig

        // new estimate for the relative wave number k*h 
        kh -= f / df;
        err = f;
    }

    // compute final result k from the dimensionless k*h 
    *pk = kh / h;
}

